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ABSTRACT 

The mechanical properties of a neutron star crust, such as breaking strain and shear 
modulus, have implications for the detection of gravitational waves from a neutron 
star as well as bursts from Soft Gamma-ray Repeaters (SGRs). These properties are 
calculated here for three different crustal compositions for a non-accreting neutron 
star that results from three different cooling histories, as well as for a pure iron crust. 
A simple shear is simulated using molecular dynamics to the crustal compositions by 
deforming the simulation box. The breaking strain and shear modulus arc found to 
be similar in the four cases, with a breaking strain of 0.1 and a shear modulus of 
~ 10^" dynecm"^ at a density of p = 10^* g cm~^ for simulations with an initially 
perfect BCC lattice. With these crustal properties and the observed properties of 
PSR J2124-3358 the predicted strain amplitude of gravitational waves for a maximally 
deformed crust is found to be greater than the observational upper limits from LIGO. 
This suggests that the neutron star crust in this case may not be maximally deformed 
or it may not have a perfect BCC lattice structure. The implications of the calculated 
crustal properties of bursts from SGRs are also explored. The mechanical properties 
found for a perfect BCC lattice structure find that crustal events alone can not be 
ruled out for triggering the energy in SGR bursts. 
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1 INTRODUCTION 

The strength of the material that comprises neutron-star 
crusts determines the maximal strength of gravitational ra- 
diation from a single deformed neutron star (Chamcl & 
Haensel, 2008). Furthermore, yielding events of the neutron 
star crust iiavc also been associated with bursts from Soft 
Gamma-ray Repeaters (SGRs) (Chamel & Haensel, 2008). 
The neutron star crust consists of the outer regions of the 
star ranging from a density of lO^gcm^'' up to nuclear den- 
sity of 2 X lO^^gcm"^ (Lattimcr & Prakash, 2004). The 
less dense regions of the crust consist of a lattice of nu- 
clei, with the ground state of a body centred cubic (BCC) 
lattice (Chamel & Haensel, 2008). At higher densities the 
crust is composed of a mixture of nuclei, electrons, protons 
and free neutrons. The mechanical properties of the neu- 
tron star crust which are important for the consideration of 
crustal deformation include the breaking strain, or maximal 
degree of deformation before yielding, the amount of stress 
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associated with this strain, as well £is the shear modulus of 
crustal material. 

Initial estimates of the breaking strain of the neutron 
star crust were made by comparing to terrestrial matter. 
Arguing that the impurities expected in the neutron star 
crust would result in weakening the structure, the breaking 
strain of the crust could range between 0^ = 10~® and = 
10~^ (Smoluchowski, 1970). The shear modulus of the crust 
lias been predicted by treating the neutron star crust as a 
Coulomb lattice, which has the approximate relationship of 
H oc (Ze)^n*/^, where the charge Z ranges between 30 and 
50, and n is the number density, which leads to a expected 
shear modulus of /i = 10'^"dyne cm^'^ (Smoluchowski, 1970). 

Recently, with molecular-dynamics simulations, the 
breaking strain of a system representing an accreted neutron 
star crust has been calculated. The crustal composition of 
the accreted crust consists of isotopes ranging from Z = 8 
to Z = 47 (Gupta et al., 2007). The breaking strain for the 
accreted crust was found to occur around 0m = 0.1, which 
is larger than had been predicted, and this result was also 
found to only moderately affected by the introduction of im- 
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purities, defects and grain boundaries (Horowitz & Kadau, 
2009). 

In this work molecular-dynamics simulations are used in 
order to investigate the crustal properties of a non-accreting 
neutron star. The molecular-dynamics calculations are car- 
ried out with the open source software LAMMPS (Large 
Atomic/Molecular Massively Parallel Simulator) (Plimpton, 
1995), which is available from Sandia Laboratories^. The 
breaking strain and shear modulus found with these simula- 
tions arc used in order to place limits on the strain amplitude 
for a gravitational wave signal expected for a fully deformed 
neutron star. Limits are also placed on the fracture size re- 
quired for the energy in a corresponding SGR burst. In the 
following section the implementation and parameters of the 
molecular-dynamics simulations are discussed. In Section 3 
the results of these simulations are presented. The simula- 
tion results are applied to the emission of gravitational waves 
and SGR bursts in section 4. In Section 5 we summarize our 
findings. 

2 COMPOSITION AND MOLECULAR 
DYNAMICS 

In order to calculate the composition of a non-accreting 
neutron star crust we used the reaction rate network 
torch^(Timmes et al., 2000). The torch software calculated 
the abundances of 489 different isotopes as the neutron star 
material cooled at densities ranging from 10® — 10^^ gcm~^. 
The compositions were calculated for three different cooling 
curves - cooling via modified Urea, cooling appropriate for 
a thick crust (Lattimer et al., 1994), as well as a thin crust. 
Each of the cooling scenarios resulted in a different crustal 
composition, see Hoffman & Heyl (2009) for more details. 

Out of the total 489 different isotopes, only the top three 
most abundant isotopes were used in the molecular dynam- 
ics simulations. For the modified Urea cooling case the top 
throe isotopes wore found to be ^"Fe, •'''*Fe, and ®°Ni. In the 
case of neutron star cooling with a thick crust the top three 
isotopes were ^''Fe, ®°Ni, and ^^Cr. In the tliiu crust case 
^''Fe, ^®Ni, and ^^Fe were found to have the top abundances. 
Table 1 summarizes these top three isotopes, as well as the 
rescaled mass fractions in order for the composition to sum 
to unity. The table also includes the number fraction of each 
isotope used in the simulations and the properties of each of 
the isotopes. 

The mechanical properties of the neutron star crustal 
material are investigated by applying a shear to the material. 
The shear is simulated using molecular dynamics. Molecular 
dynamics is the integration of Newton's equations of motion 
for a system of particles. The molecular dynamics simula- 
tions calculate the effect of a specified force applied collec- 
tively to a system. Within the simulations the position, ve- 
locity, and acceleration of each particle are calculated, lead- 
ing to a determination of the thermodynamic properties of 
the system. For the systems of interest in this work the par- 
ticles interact via a Yukawa, or screened Coulomb, potential 

ZiZkC —K.r /IN 

Vij = — - — e , (1) 
^ http://lammps.sandia.gov 
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Isotope 


A 


Z Mass Rescaled 


Number 


Mass 






Fraction X 


Fraction 




Modified Urea Composition 


56 pg 


56 


26 0.559 0.6014 


0.5961 


1.0 


54Fe 


54 


26 0.3187 0.3429 


0.3524 


0.964 




60 


28 0.05175 0.0557 


0.0515 


1.071 


Thick Crust Composition 


56 


.56 


26 0.9286 0.9429 


0.9431 


1.0 


60Ni 


60 


28 0.03138 0.0319 


0.0298 


1.071 


B2Cr 


52 


24 0.02485 0.0252 


0.0271 


0.929 


Thin Crust Composition 


54Fe 


54 


26 0.6477 0.6675 


0.6808 


0.964 


58Ni 


58 


28 0.2209 0.2276 


0.2161 


1.036 


56 Pg 


56 


26 0.1018 0.1049 


0.1031 


1.0 



Table 1. The top three isotopes from the torch calculations for 

three difEcrcnt cooling scenarios: modified Urea, thick crust, and 
thin crust cooling. The original mass fraction from the torch cal- 
culations, as well as the rescaled mass fraction (X) are included. 
The number fractions determine the number of each type of iso- 
tope used in the molecular dynamics simulations. The mass value 
is the scaled simulation parameter, where in the simulations the 
mass is unity for A = 56. 

where r is the distance between the particle pair, Z is the 
particle charge, and k is the inverse screening length. The in- 
verse screening length is given as k = 2a^^'^kF /'k^^'^ , where 
kp ~ (STT'^ne)^''''*, with He — {Z)n as the electron num- 
ber density. The molecular dynamics sirrmlations are con- 
ducted using LAMMPS, with the evolution of the particle 
motion calculated with velocity- Verlet integration (Swope 
et al., 1982). The velocity- Verlet algorithm is a variation of 
Verlet integration (Verlet, 1967). It is a second-order inte- 
gration method which is time reversible and is a result of 
the expansion of the particle position with time. At a given 
timestep the velocity- Verlet algorithm calculates both the 
velocity and position of each particle, as well as the force of 
each of the interactions. In order to decrease the calculation 
time, LAMMPS makes use of parallel computing (Plimpton, 
1995) and neighbour lists. The neighbour lists each contain 
a list of particles within a specified distance from each other 
and only particle pairs within a specified cut-off radius are 
used in the interaction calculations. LAMMPS version 10 
Jul 2010 was used in the simulations discussed here. 

The simulation parameters are expressed in reduced 
units based on the Yukawa potential (Eq. 1). These reduced 
units introduce a set of characteristic quantities: length, en- 
ergy, mass, and charge. The characteristic length scale, a is 
a = n~^^^, where n is the number density. The characteristic 
energy scale is defined as 



which is expressed with respect to a charge of Z = 26. The 
characteristic mass, nia, is expressed in units of the mass 
of iron, or rria ~ A/ 56, thus for an iron ion ma is unity. 
In a similar manner, the characteristic charge is expressed 
in units of the iron charge, thus Za = Z/26. These charac- 
teristic quantities are used in order to relate reduced unit 
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Unit Reduced Unit 



Temperature 
Energy 


E* 




Distance 


r* 




Pressure 


P* 


= p — 

Ua 

-(^) 


Time 


f -- 



Table 2. The reduced unit system, based on the Yukawa po- 
tential, used in tiic LAMMPS simulations. Parameter values quoted 
from the simulations are in reduced units. The physical quanti- 
ties are determined by using appropriate values the characteristic 
distance, a, and energy, Ua, parameters for the material. 





56Fe -'^Cr 


54 




60Ni 


sepe 


1.0 0.9231 


1.0 


1.0769 


1.0769 


52Cr 


0.8521 


0.9231 


0.9941 


0.9941 


54Fe 




1.0 


1.0769 


1.0769 








1.1598 


1.1598 


60Ni 








1.1898 



Table 3. The pair-coefficients, ZjZ, e^, for the different parti- 
cle pairs. The pair-coefficients are used in the calculation of the 
interaction potential. These pair-coefficients are used in the four 
different cases of crustal simulations, the three impure cases and 
a pure iron crust case. 







(Z) 


K 


rc 


Qimp 


Mod. Urea 

Thick 
Thin 


55.501 
56.011 
55.071 


26.103 
26.005 
26.432 


0.8847 
0.8836 
0.8884 


9.043 
9.054 
9.005 


0.195 
0.228 
0.678 



simulation parameters to the corresponding physical quan- 
tities. For example the characteristic quantities of ^^Fe at 
a density of p = lO^^g cm~^ are a characteristic length of 
a ~ 10~^^cm and a characteristic energy of Ua ~ 10~^erg. 
The relation between the reduced and physical values are 
listed in Table 2. 

Within all the simulations the the number density is set 
to unity. Periodic boundary conditions are also used with the 
shear introduced by deforming the box in the x-direction. 
The degree of strain is determined by taking the ratio of the 
displacement of the top of the simulation box to the original 
box length, Ax/l. The strain rate is the the velocity of the 
shear applied at the top of the box divided by the length of 
the box, v/l. 

The simulations all begin with each of the particles 
given an initial velocity kick, which gives the temperature 
of the simulation below the melting temperature. After 200 
steps of equilibrating the particles to a specified tempera- 
ture the box is deformed, holding the temperature constant 
through temperature rescaling. The velocity- Verlet integra- 
tion, which calculates the positions and particle velocities, is 
performed on a system of particles which is consistent with a 
microcanonical ensemble, constant number, volume and en- 
ergy (NVE). A nearly constant temperature is maintained 
by rescaling the velocities. A system of particles consistent 
with a canonical ensemble, constant number, volume, and 
temperature (NVT), could have also been selected and used 
for temperature control, but the velocity-rescaling thermo- 
stat controls the temperature in a more gentle and direct 
manner than the Nose-Hoover or Langevin thermostat. Of 
course, the latter would better approximate the canonical 
ensemble (NVT); however, for simplicity and for comparison 
with previous work (e.g. Horowitz & Kadau, 2009), we also 
used velocity rescaling. The interaction calculation is cut-off 
at a distance of rc — S/n between particle pairs, where k 
is the inverse screening length. The pair-coefficients of the 
different particle pair interactions, ZiZjC^, arc listed in Ta- 
ble 3. The average mass, average charge, inverse screening 
length, cut-off radius, and impurity factor of the composi- 
tion are listed in Table 4 for the three impure cases. 



Table 4. The inverse screening length and the corresponding cut- 
off radius used in the simulations. The parameters {A) and {Z) 

were used in the calculation simulation parameters. The Qimp 
value, where Qi^p = (2^) - (Z)'^ (Itoh & Kohyama, 1993), indi- 
cates the degree of impurity of the sample for the three impure 
crustal cases. 



3 RESULTS 

The initial molecular dynamics simulations performed in- 
cluded tests to examine the effect of the size of the simula- 
tions as well as varying the strain rate in order to look for 
convergence of the material properties examined. In order 
for the shearing simulations to be run below the melting tem- 
perature, the melting temperature of the various materials 
was also determined. With these initial test simulations, pa- 
rameters for the simple shear simulations were determined. 
The simple shear simulations were all performed at a tem- 
perature below the melting temperature of the crystal. A 
simulation box size of 25 x 25 x 25 unit cells for simple shear 
simulations of 31250 particles. Strain rates of 20x, 5x, and 
2.5x10"'' were examined in the simple shear calculations. 
The application of the shear was also used on crystals with 
an imperfect crystal structure and compared to the results 
of a perfect BCC crystal lattice structure. The temperature 
dependence of the stress-strain relationship was also exam- 
ined for three different temperature regimes. The possibility 
of a second yielding event was investigated via simulations 
run to a larger degree of strain. Finally, the reversibility of 
the yielding event, due to the applied shear, was also exam- 
ined. 



3.1 Test Simulations 

The test simulations included determining the melting tem- 
perature of the different compositions as well as investigat- 
ing size effects. The melting temperature was determined in 
order to ensure the simple shear simulations were below the 
melting point of the compositions. The box size tests were 
performed in order to determine the effect the box size had 
on the simulation, as well as finding the smallest number of 
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particles which could be used to gain a reliable result from 
the simulations. 

3.1.1 Melting Temperature 

The melting temperature of the three impure crustal compo- 
sitions, as well as the pure iron case, was determined by in- 
creasing the temperature of the simulation box. In these sim- 
ulations, unlike the shearing simulations, the temperature 
was controlled by a Langevin thermostat. With a Langevin 
thermostat the system of particles is within a heat bath 
and the temperature is controlled through a frictional force 
as well as a random force, which imparts a random veloc- 
ity kick (Schneider & Stoll, 1978). The melting temperature 
was defined as the temperature in which the mean squared 
displacement of the particles has the greatest change in the 
value between time steps. The simulation boxes were heated 
from a temperature of T* = 0.001 to a value of T* = 0.02 
over 10** steps. All four cases were found to melt at around 
the same temperature of T*~0.01. To compare this to a 
physical temperature, at a density of p = lO^"' g/cm"* these 
melting temperatures correspond to ~ 7 x 10** K. The tem- 
perature of the simple shear simulations are always less than 
this melting temperature. 

3.1.2 Size Effects 

In order to examine the possible simulation size effects a 
BCC iron lattice with periodic boundary conditions was 
strained at a rate of 20 x 10~^ for five different box sizes. The 
simulation box sizes included boxes with unit cell volumes 
of 10 X 10 X 10, 12 X 12 X 12, 20 x 20 x 20, 25 x 25 x 25, and 
35 X 35 X 35, corresponding to 2000, 3456, 16000, 31250, and 
85750 particles, respectively. The result of deforming these 
five different simulation box sizes are compared in Figure 1, 
which displays an apparent difference in the box size and 
yielding stress and strain. The shear modulus is the same in 
all five cases, but simulations with a larger number of parti- 
cles are found to yield earlier. The stress-strain relationship 
for the 35 X 35 X 35 and 25 x 25 x 25 unit class box sizes 
are similar, and in order to decrease the calculation time a 
25 X 25 X 25 box size is used for simulations of the different 
crustal compositions. 

3.2 Simple Shear 

The four cases of three impure compositions appropriate for 

cooling by modified Urea, a thick crust, and a thin crust, 
as well as a pure iron crust were deformed in molecular dy- 
namics simulations in order to determine the mechanical 
properties. The shear was applied by deforming the simu- 
lation box in the x-direction. Three different strain strain 
rates were applied in the simulations, 20, 5, and 2.5x10"^, 
in order to test for convergence of the mechanical proper- 
ties. The simulation boxes had periodic boundary conditions 
and contained 31250 particles in all cases. Each simulation 
stared in a perfect BCC lattice configuration, for the im- 
pure cases the isotopes were randomly assigned to lattice 
locations. The stress-strain relationships at a stain rate of 
20 X 10~® for the four cases are compared in Figure 2. The 
yielding strain and shear modulus determined for the four 
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Strain 

Figure 1. The stress-strain relationship of a pure BCC crystal 
comparing five different sizes of simulations. The legend indicates 
the unit cell volume of the simulation box. The number of par- 
ticles range from 2000 in the smallest simulation box to 85750 
particles in the largest simulation box size. For the same applied 
strain rate of s = 20 X 10~® the five sizes of simulations share 
the same shear modulus, but the breaking strain is different. The 
larger simulation boxes are found to yield before the smaller sim- 
ulation boxes. 

cases at the three different strain rates are listed in Table 
5. The pure iron, modified Urea and thick crust composi- 
tions share very similar breaking strains, of (f>m ~ 0.11, and 
shear moduli, of n* ~ 0.21. The thin crust composition is 
found to yield to the applied stress at a later degree of de- 
formation than the other three compositions, at around a 
breaking strain of 4>m ~ 0.12. These results of a breaking 
strain of ~0.1 for non-accreting crustal compositions is ten 
times larger than the initial estimates from Smoluchowski 
(1970), but the results arc consistent with previous molecu- 
lar dynamics simulations for an accreted crustal composition 
(see Horowitz &: Kadau, 2009). 

3.3 Imperfect Crystal Structure 

The shearing simulations discussed above initially had a per- 
fect BCC lattice structure. The addition of defects to the 
simulation lattice would be expected to decrease the crystal 
strength, as a result the perfect BCC lattice structure simu- 
lations may represent an upper limit on the crustal strength. 
It should be noted that previous molecular dynamics simula- 
tions which included defects were found to only moderately 
affect the strength (Horowitz & Kadau, 2009). In the case of 
a non-accreting neutron star crustal composition the effect 
of defects are examined by melting the perfect BCC lattice 
structure and then allowing the simulation box to recrystal- 
lize. The application of a shear to a perfect BCC lattice and 
an imperfect lattice is compared in Figure 3 for a thick crust 
composition. The behaviour, as displayed in Figure 3, of the 
two different initial crystal types have different behaviours 
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Figure 2. The stress-strain relationship for tiic tiucc impure 
crustal compositions: modified Urea (mu3), thick crust (thick), 
and thin crust (thin), compared to a pure iron BCC lattice (bc- 
cfe). A strain rate of 20 x 10~® was used in all four of these 
simulations. The shear modulus was found to be similar for all 
four of the simulations. The thin crust yielded to the deforma- 
tion at a higher strain than the other three simulations. As these 
results use a perfect BCC lattice structure, these simulations set 
an upper limit for the shear modulus and breaking strain of a 
neutron star crust. 



Composition 


Strain Rate 


Yielding 


Shear 




(xlO-6) 


Strain 


Modulus 


Pure 


20 


0.114508 


0.27217 


Iron 


5 


0.112896 


0.27219 




2.5 


0.111295 


0.27220 


Modified 


20 


0.115620 


0.2735 


Urea 


5 


0.113396 


0.2736 




2.5 


0.113174 


0.2735 


Thick 


20 


0.113285 


0.27179 


Crust 


5 


0.111895 


0.27178 




2.5 


0.111673 


0.27179 


Thin 


20 


0.122846 


0.26644 


Crust 


5 


0.120122 


0.27620 




2.5 


0.119455 


0.27711 



Table 5. A comparison of the the yielding strains and shear mod- 
uli for the four different compositions. The ratio of the strain 
rate to the plasma frequency is 56.42, 3.979, and 1.989x10"^, for 
strain rates of 20, 5, and 2.5xl0~®, respectively. 



Figure 3. A comparison of the stress-strain relationships of a per- 
fect and an imperfect crystal which contains defects for the case 
of a thick crust composition. A strain rate of s = 20 x 10~® was 
applied in the X-direction for both simulations. The two crystals 
display differing behaviours in their stress-strain relationships. 



3.4 Temperature Dependence 

The temperature dependence of the stress-strain relation- 
ship was investigated by performing additional simulations 
at temperatures twice and half as hot as the original sim- 
ulation, for a total of three different temperatures exam- 
ined. Thus, for each crustal composition, temperatures of 
T* = 0.002, T* = 0.001, and T* = 0.0005, which corre- 
spond to physical temperatures of 6.9 x 10* K, 3.5 x 10® K, 
and 1.7 x 10* K at a density of p = lO^** gcm~^, were inves- 
tigated. A strain rate of 20 x 20^'' was applied to an initially 
perfect BCC lattice structure in all cases. The three temper- 
ature stress-strain relationship results for a pure iron perfect 
BCC lattice are compared in Figure 4. Though the rela- 
tionships are very similar, the shear modulus for the colder 
simulations is slightly larger than the hotter temperatures 
and the material yields to applied stress earlier in the colder 
simulations as opposed to the hotter. This behaviour is also 
observed in the three impure crustal compositions, with the 
shear modulus and the breaking strains listed in Table 6 for 
the four cases investigated. The temperature dependence ob- 
served in these simulations is expected from models of shear 
modulus-temperature relationships, e.g. Nadal & Le Poac 
(2003), where the shear modulus is found to decrease with 
increasing temperature. 



in terms of the stress-strain relationship and the resulting 
shear modulus and the yielding strain. The difference in the 
effect of defects as reported here as opposed to Horowitz & 
Kadau (2009) could be due to the number of added defects 
or grain boundaries. 



3.5 Second Yielding Events 

In the lifetime of a neutron star, the crust would have 
likely undergone multiple yielding events. Subsequent yield- 
ing events, after the first major event, are investigated in 
these simulations by deforming the simulation box to a 
strain of 0.4, where the original simulations were deformed 
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Figure 4. The stress-strain relationship for pure iron BCC crys- 
tals at three different temperatures. At T*~ 0.01 the crystal 
melts; all simulation temperatures are below the melting temper- 
ature. A strain rate of 20 X 10^® was applied in all three cases. For 
each temperature the shear modulus and breaking strain slightly 
differ, though T* = 0.001 and T* = 0.0005 were found to be 
very similar. The shear modulus of the colder simulation was 
slightly larger then the warmer simulations. The colder simula- 
tion is found to yield at a smaller strain, but higher stress, than 
the hotter simulations. 



Crust 


Temperature 


Breaking 


Shear 


Composition 


T* 


Strain 


Modulus 


Pure 


0.002 


0.11712 


0.26069 


Iron 


0.001 


0.114508 


0.27217 




0.0005 


0.112285 


0.27821 


Modified 


0.002 


0.116064 


0.26253 


Urea 


0.001 


0.115620 


0.27351 




0.0005 


0.112952 


0.29558 


Thick 


0.002 


0.115397 


0.26052 


Crust 


0.001 


0.113285 


0.27179 




0.0005 


0.112396 


0.27789 


Thin 


0.002 


0.122846 


0.26644 


Crust 


0.001 


0.121901 


0.27637 




0.0005 


0.122401 


0.28060 



Table 6. The temperature, breaking strain and shear modulus 
for the three different impure crustal compositions, as well as the 
pure iron case. The three different crust compositions melt at 
approximately T* = 0.01, all temperature simulations are below 
the melting temperature. At a density of p = lO^^g/cm^ the 
corresponding physical temperatures are 6.9 x 10* K for T* = 
0.002, 3.5x10* KforT* = 0.001 and 1.7 x 10* K for T* = 0.0005. 



Figure 5. The stress-strain relationship of a 31250 particle, pure 
iron BCC crystal with an inverse screening length of ft = 0.8835 
with a strain rate of 20 x 10^® applied. The simulation box was 
deformed to a strain of Ax/i = 0.4 in order to investigate the 
occurrence of a second material failure. A second yielding event 
occurs for the crystal between a strain of 0.133 to 0.186 with a 
shear modulus of ;U*~ 0.11 in reduced pressure units. 



to a 0.2 stain. In these longer simulations a strain rate of 
20 X 10^® was applied to all four types of crustal compo- 
sition. The simulations discussed here started with a per- 
fect BCC lattice structure. The stress-strain relationship for 
the pure iron case is displayed in Figure 5. As the stress is 
applied, the crystal undergoes a linear regime before yield- 
ing to the applied stress, but after the major yielding event 
the stress-strain relationship undergoes an additional linear 
regime, which has a different slope, or shear modulus, as the 
initial event. In the case of the pure iron BCC crystal the 
shear modulus before the initial yielding point is /i* ~ 0.27 
and in the second linear regime, measured between values of 
a strain of 0.133 and 0.186 the shear modulus is ~0.11. 
In the cases with the modified Urea, thick and thin crustal 
compositions similar behaviour occurs with the shear mod- 
ulus changing after the initial yielding point. In the case 
of the thick crust composition the reduced shear modulus 
changes from /x* ~0.27 before the first major yielding event 
to [i* ~ 0.12 for the second. Comparing to the results of the 
imperfect crystal in section 3.3, the behaviour of the second 
break is similar to the imperfect crystal, though the linear 
regime has a slightly larger slope for the initially perfect 
lattice structure crystal. 



3.6 Yielding Reversibility 

In order to investigate the reversibility of the yielding event, 
the direction of the box deformation was reversed after the 
major yielding event occurs. If the yielding event was re- 
versible than the shear modulus would be the same deform- 



Mechanical Properties 7 



0) 

an 
0) 



_l 1 1 1 1 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 i_ 


: / 


Reverse 

- 


/ 


- 











'0.00 



0.05 



0.10 0.15 

Strain 



0.20 



Figure 6. The stress-strain relationship of a pure iron BCC crys- 
tal with an inverse screening length of k = 0.8835 and a strain 
rate of 20 x 10~® applied. The deformation of the box was reversed 
after the material yields. For illustrative purposes the solid Une 
indicates the strain as is the box deformation was continued in the 
forward direction, the dotted line indicates the box deformation 
as it is being deformed in the reverse direction. In the forward 
direction the shear modulus was fi*/^ 0.27, after the box defor- 
mation was reversed the shear modulus is 0.17. The shear 
modulus was measured between a strain of 0.175 and 0.2115. 



ing in the reverse direction as the the forward direction. In 
these simulations the simulation box containing a perfect 
BCC lattice is deformed at a strain rate of 20 x 10^^, upon 
yielding the box is then deformed in the opposite direction, 
but at the same strain rate. The stress-strain relationship 
for this type of simulation is displayed in Figure 6 for a pure 
iron perfect BCC lattice structure; for plotting purposes the 
amount of strain in the figure continues to increase after re- 
versing the shear direction as (step sizexAt x 20 x 10~®). 
In this case of the pure Iron BCC crustal the shear modu- 
lus is fj,* ~0.27 in the forward direction. After reversing the 
direction of the shear direction the shear modulus is mea- 
sured to be ^*~0.17, measured between a strain of 0.175 
and 0.2115 on the corresponding figure. As the shear modu- 
lus does not return to the original value, the yielding event 
is not reversible, similar results are found for the three im- 
pure composition cases. The reduced shear modulus in the 
reverse direction, /i* ~0.17, is larger than what was found 
for the second break, /j* ~ 0.11 in the pure iron crystal case. 
Consequently, the shear modulus for the reverse direction 
would also be expected to be greater than for the imperfect 
crystal. 



4 DISCUSSION 

The simple shear simulations of both the non-accreted 
crustal compositions, presented above, and the accreted 
crustal composition (see Horowitz & Kadau, 2009) find that 



the neutron star crust yields at a strain of ~0.1, indicating 
that the neutron star crust is much stronger than the initial 
estimates. The results of a stronger crust could have implica- 
tions for the observational consequences of crustal yielding 
events. Comparing the results of the pure iron case to the 
simulation results which contained three different isotopes 
did not indicate much of a change in the material properties 
for simulations with impurities and those with out. The ad- 
dition of imperfections, while the system may globally have 
a BCC structure, indicate a difference in the stress-strain 
relationship as compared to systems which started in a per- 
fect BCC lattice structure. Depending on the strength of the 
neutron star crust, detectable gravitational waves may be 
emitted due to mountains supported on the surface. Yield- 
ing events in the crust have been associated observationally 
with bursts from SGRs. The results of the shearing simu- 
lations on a perfect BCC lattice structure are now applied 
to these two observational consequences of the mechanical 
properties of the neutron star crust. 

The strain at which the crustal material was found to 
yields was at ~0.1, and the shear modulus from the simula- 
tions was found to be 0.27 in reduced pressure units. To 
apply these results to gravitational waves and SGR bursts 
the simulation parameters need to be converted to physi- 
cal units. The reduced pressure, or the pressure from the 
simulations, is expressed as: 

„3 



P* = P 



Ua 



(3) 



where P is the physical pressure, Ua is the characteristic en- 
ergy, and a is the characteristic distance. The observational 
consequences of crustal cracking are typically investigated 
at the bottom of the neutron star crust, at a density of 
p = 10^*gcm~^. These simulation results are only appro- 
priate for densities closer to the neutron star surface, but 
in order to compare to previous work, the results from the 
simulations are extrapolated to the greater density by tak- 
ing into account the neutron fraction. At the bottom of the 
crust the neutron fraction is Xn = 0.8, which is higher than 
at lower densities closer to the neutron star surface. The 
charge can be set to Z = 20 and the atomic mass to ^ = 88 
at the bottom of the crust (Ushomirsky et al., 2000). With 
these considerations the number density is related to the 
neutron fraction as: 



n = a-^ = ^(l 
m 



X„) = 



10"g/cm3 



X 0.2, 



(4) 



10^* g cm ^, the bottom of the neu- 



88 rriu 

where niu is the atomic mass unit. 

At a density of p ■ 
tron star crust, the characteristic length is a = n~^^^ = 
1.94 X 10~^'^cm or 19.4 fm and the characteristic energy, 
where the charge is Z = 20, is Ua = {Zef/a = 4.76 x 
10~^ erg, which corresponds to ~30 MeV. With these char- 
acteristic values, the shear modulus of /u* = 0.27 from the 
simulations corresponds to p, — 1.76 x 10^" dyne cm~^ at 
the bottom of the crust. This value is close to the perviously 
used value of 10^" dynecm^'^ (for example Ruderman, 1969). 
Note that at a density of p = lO'^ g cm~^ the shear modulus 
would be /X ~ 10^^ dynecm"^. The strain at which yielding 
occurs is a unitless quantity, thus the simulation values are 
the physical values with a breaking strain of 4>m ~ 0.1 used 
in order to compare the simulation results to observations. 
Note that the breaking strain found for the non-accreted 
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neutron star crust is in agreement of the breaking strain cal- 
culated in Horowitz & Kadau (2009) of close to 4>m = 0.1, 
for an accreted crust crustal composition. 



4.1 Gravitational Waves 

Gravitational waves may be emitted by a neutron star with 
varying quadrupole moment. This quadrupole moment is 
the result of the neutron star perturbed to be asymmetric 
about its rotation axis due to a mountain, or deformation, 
on the neutron star surface. The deformation the crust can 
sustain depends on its mechanical properties, such as the 
breaking strain. The maximum quadrupole can be expressed 
in terms of the breaking strain and shear modulus of the 
crust (Ushomirsky et al., 2000). In these calculations the 
breaking strain is expressed as Q22 = 7 4>m I, where 7 is a 
factor which depends on the mass and radius of the neutron 
star, 0m is the breaking strain, and I is an integral which 
is directly proportional to the shear modulus, and thus Q22 
is directly proportional to the shear modulus, jj,. Scaling by 
the value of Q22 for a shear modulus of /x = 10^° dyne/cm^ 
yields the following relationship 

^°"^-^(lF^)( lO3Ody:ecm-0 - 

For the simulation results of the perfect crystal structure 
with a breaking strain of 0m ~ 0.1 and at a density of 
p = 10^^gcm~^ where = 1.76 x lO^^dynecm'^ the 
quadrupole is Q22 = 1.76 x lO^^gcm^, assuming the crust 
is maximally deformed. 

The effect of a mountain on a neutron star is observed 
through the measured strain amplitude of gravitational ra- 
diation, which is given by 

, 16 /7r\i/3 GQ22 

where is the angular frequency and d is the distance to 
the neutron star (Ushomirsky et al., 2000). Observational 
upper limits have been placed on the gravitational wave 
strain amplitudes of 78 radio pulsars using the Laser Inter- 
ferometer Gravitational Wave observatory (LIGO) (Abbott 
et al., 2007). A detection of gravitational waves from a neu- 
tron star, with a comparison to predicted strain amplitudes 
could possibly put constraints on the strength and structure 
of the neutron star crust. 

Of the 78 different pulsars which have had upper lim- 
its placed on their gravitational wave emission, the prop- 
erties of three pulsars are examined in the context of the 
simulation results: PSR J1603-7202, PSR J2124-3358, and 
PSR J0534-F2200 (the Crab pulsar). Though these pulsars 
have likely undergone accretion events as two of the pul- 
sars are recycled and the third is the Crab pulsar, which 
is in a supernova remanent and may have accreted some 
of the nearby matter, the breaking strain results found for 
the non-accreted crust are similar to those reported for an 
accreted crust in Horowitz & Kadau (2009). The LIGO ob- 
servational limits of the gravitational wave strain amplitude, 
h, are compared to the strain amplitude for the deformation 
of an initially perfect BCC lattice crust which is maximally 
deformed in Table 7. The distance, d, and the spin frequency, 
v, of the three pulsars are also listed in the table. 



Pulsar 


d 


1/ 


LIGO 


Predicted 




kpc 


Hz 


Max 


Max 








log(/i) 


log(ft) 


PSR J1603-7202 


1.64 


67.38 


-24.58 


-25.8 


PSR J2124-3358 


0.32 


202.71 


-24.31 


-24.1 


PSR J0534+2200 


2 


29.80 


-23.51 


-26.6 



Table 7. The distance, d, spin frequency, u, and the obser- 
vational and predicted gravitational wave strain amplitude, h, 
for the three pulsars of interest. The predicted gravitational 
wave strain amplitudes are below the observational LIGO up- 
per limits for two of the cases. For pulsar PSR J2124-3358 the 
predicted strain amplitudes indicate that if the crust was ini- 
tially perfectly BCC in structure and maximally deformed, than 
gravitational waves should have been detected from this source. 
The distances of each pulsar were taken from the ATNF Pulsar 
Catalogue (http : //www. atnf . csiro . au/people/pulsar/psrcat 
(Manchester et al., 2005)). The spin frequency and LIGO up- 
per limits are from Abbott et al. (2007) . The predicted maximum 
strain amplitudes are calculated using the distance and spin fre- 
quency of each, of the pulsars, as well as the quadrupole as calcu- 
lated using the results from the simulations. 

With the mechanical properties calculated in the molec- 
ular dynamics simulations it is found that if PSR J2 124-3358 
was maximally deformed than it would have been detected 
as a source of gravitational waves. Assuming that the crustal 
strain is uniforndy distributed to the breaking point, the 
non-detection of gravitational waves from this neutron star 
suggests that either the neutron star crust is not a perfect 
BCC crystal, or that the crust is not maximally deformed. 
The search for gravitational waves from pulsars could be 
used as a method to determining the structure, strength, or 
degree of deformation of the neutron star crust. 

4.2 Bursts from Soft Gamma-ray Repeaters 

Soft Gamma-ray Repeaters (SGRs) are characterized by 

their short bursts of low-energy gamma rays, which can 
reach a peak luminosity of 10'*^ erg s^^ (Woods & Thomp- 
son, 2006). There have been eleven SGRs observed; seven 
confirmed and four candidates^. Of these eleven SGRs three 
have been also observed to exhibit giant flare events, which 
can reach luminosities of 10'*'* erg/s (Hclfand & Long, 1979). 
The observed SGR bursts have been associated with crustal 
activity; the bursts may be triggered by the yielding events 
due to stresses placed on the crust by the magnetic field 
(Chamel & Haensel, 2008). The strength of the magnetic 
field required to trigger a burst is dependent on the break- 
ing strain of the crust: 

where Esgr is the energy of the burst, / is the length of the 
crustal fracture and 4>m is the breaking strain (Thompson & 
Duncan, 1995). Using the results from the simulations, the 
observed energy emitted during burst and a neutron star's 
magnetic field strength, limits can be placed on the size of 

^ Data from the McGill SGR/AXP online catalog, 
http : //www . physics . mcgill . ca/- pulsar /magnetar /main . html 
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the crustal fracture. By estimating the fracture size, this can 
be used as an indication of if the burst can be attributed 
to the crust alone. If the fracture length required for an 
observed energy of a burst is larger than the crust size, then 
the burst can not be due to crustal events alone. 

The length of the fracture, using Equation 7, can be 
expressed as: 

For a typical burst energy of Esgr = W'^^ erg, assuming 
a magnetic field of S = 10^^ G, and using the simulation 
result of a yielding strain of 0.1 the corresponding fracture 
size is Z = 10 m. In the case of the more energetic giant flares 
with Esgr = 10** erg the corresponding fracture length is 
I = 300 m. The neutron stax crust has a depth of approx- 
imately 1km with the outer crust encompassing ~100m. 
The fracture length associated with typical burst energies, 
/ = 10 m, is within the confines of the crust, thus these events 
could be attributed to crustal activity. The giant flare events 
require a larger fracture length, / = 300 m, which would en- 
compass ~15 — 30% of the crustal depth. This larger fracture 
could possibly propagate along the surface or vertically and 
is within the size constraints of the neutron star crust. 

The crustal fracture lengths were calculated using prop- 
erties of a perfect BCC lattice structure. The simulations 
with an imperfect crystal structure were found to dis- 
play different mechanical projjcrties from the perfect crystal 
structure simulations. A neutron star crust which yields at 
smaller strains would require a larger fracture length. De- 
ponding on the structure of the neutron star crust, how close 
it is to a perfect BCC lattice, crustal events could be ruled 
out as the mechanism for SGR bursts. 



5 CONCLUSIONS 

In order to investigate the mechanical properties of a 
non-accreting neutron star crust we have performed var- 
ious molecular dynamics simulations using the software 
LAMMPS. With these simulations tests were first performed 
in order to test the size effects of the simulation box size, 
as well as determining the melting temperature of the var- 
ious crustal compositions. A simple shear was applied to 
the crustal compositions by deforming the initially perfect 
BCC lattice simulation box of four different crustal compo- 
sitions, those appropriate for a neutron star cooled via mod- 
ified Urea, with a thick crust, and with a thin crust, as well 
as a pure iron crust. Additional simulations were performed 
to compare the stress-strain relationship of a system initially 
with a BCC lattice structure to that which included defects. 
The temperature dependence on the the breaking strain as 
well as the shear modulus was also investigated with ad- 
ditional simple shear simulations at different temperatures. 
The possibility of a second yielding event as well as the re- 
versibility of the yielding event were also investigated in the 
four cases. The calculation of the breaking strain and shear 
modulus from the simulations was applied to the context of 
gravitational wave emission from a neutron star due to a 
deformation on the neutron star surface, such a mountain. 

The mechanical properties for the different crustal com- 
positions, which include pure and impure compositions, with 



an initial BCC lattice structure were found to share similar 
characteristics, this includes the pure and impure simula- 
tions. The shear modulus at a density of p = 10^* g cm~^ was 
found to be around 1.78 x 10^° dyne cm~^ and the breaking 
strain was found to be ~ 0.1 in all four cases. This breaking 
strain measurement is in agreement with that found for an 
accreted neutron star crust composition (Horowitz & Kadau, 
2009). As the shear was only applied in one crystallographic 
orientation, deforming along a different axis could result in 
different measured properties. With these mechanical prop- 
erties the constraints on the detection of gravitational waves 
were considered. The properties of three pulsars were con- 
sidered in order to predict the strain amplitude associated 
with a neutron star with the calculated mechanical proper- 
ties. LIGO had placed upper limits on the gravitational wave 
emission of 78 different pulsars and of the three investigated, 
one was found to have a strain amplitude larger then the up- 
per limit placed on the neutron star. This indicates that if 
the crust had these calculated properties and was maximally 
deformed, then gravitational waves would have already been 
observed with LIGO from PSR J2124-3358. In the case of 
bursts from SGRs the fracture lengths required for observed 
burst energies were found to be within the confines of the 
neutron star crust. 

Crustal material which included defects as compared to 

that which had an initially perfect BCC lattice structure 
were found to have differing stress-strain relationships. The 
defects in the simulations were introduced by first melting 
the perfect BCC lattice and then cooling the material to a 
solid. This behaviour of the initially imperfect crystal struc- 
ture was found to display a similar shear modulus to that of 
the perfect BCC lattice crystal after the first major yielding 
event occurs. The effect of defects on an accreted neutron 
star crust were examined in Horowitz & Kadau (2009) by 
introducing six grains to the simulation. The addition of 
defects in the accreted crust case were found to affect the 
stress-strain relationship only moderately. The differences 
between the defects in the non-accreted and accreted crust 
composition could be due to the number of grains in each 
of the simulations. Future work would include investigating 
the affect of grain rmmber on the corresponding stress-strain 
relationship. It should be noted that the crustal composition 
in the accreted case has an impurity factor which is higher 
than the limits placed on the factor through neutron star 
cooling (Horowitz & Kadau, 2009). A neutron star which 
has not undergone any accretion events is unlikely, and in 
the accreted crust case the impurity factor is higher than ob- 
servational constrains, the results from these two cases may 
represent two extremes of what may be occurring in nature. 
In a low mass X-ray binary systems with an accretion rate 
of lO^^Mg yr^^ the crust could be replaced via accretion 
in lO'^ years (Chamel & Haensel, 2008). In the case of dim 
isolated neutron stars, as well as AXPs and SGRs the accre- 
tion rate has been estimated to ranee from 3.2x 10"gs-i to 
4.2 X 10^ '^gs^^ or in solar mass units the accretion rates are 
5 X lO^^^Moyr-i to 7 x IO'^Mq yr'^ (Alpar, 2001). For 
an isolated neutron star these accretion rates correspond to 
a total crustal mass of O.OIMq (Chamel & Haensel, 2008) 
being replaced in ~10'' to ~10^years. The accretion rate 
can also be estimated by attributing the observed luminos- 
ity to accretion. For at 10 km and 1.4M0dim isolated neu- 
tron star, with a luminosity of lO^^erg, this corresponds to 
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M ~ 1O~^'*M0 yr~^ and the crust being replaced by accre- 
tion in 10^^ years. As such, the crustal composition calcu- 
lated for the non-accreting neutron stare, presented in this 
paper, would hold for those isolated neutron stars which 
have an accretion rate on the lower end of the rates. 

The molecular dynamics simulations reported here leave 
many venues for future work in order to fully understand 
the neutron star crust. The neutron star crust would not be 
expected to have a perfect lattice structure, thus it would be 
important to understand the effect of defects added to the 
simulation, as well how the orientation of the applied shear 
also affects the material. Fully characterizing the material as 
the simulation progresses would give an indication of how 
the structure of the crustal material changes with applied 
shear and cooling, this may have an effect on the material 
properties, as well as conductivity, of the crust. In all the 
simulations magnetic fields were not directly placed within 
the simulation, the addition of magnetic fields may have an 
effect on the mechanical properties, especially as with high 
magnetic fields there is a change in the electron screening 
(Sharma & Reddy, 2011). The simulation results reported 
here indicate upper limits on the crustal properties, further 
simulations may bring us closer to understanding the true 
nature of the neutron star crust. 
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